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ABSTRACT 


The  linear *  discrete  Kalman  filter  was  analysed  us ins  a 
f reouency-domain  approach.  Process  and  measurement  noise 
covariances  are  shown  to  be  critical  design  parameters  which* 
together  with  the  assumed  prior  state  and  covariance 
estimates*  completely  determine  the  Sain  schedule  of  the 
linear  Kalman  filter.  Several  relevant  desisn  techniques  are 
illustrated  and  discussed.  The  concepts  of  smoothinS  and 
sharpenins  are  demonstrated.  Extensions  to  adaptive* 
non-linear*  and  nan-parametric  filterins  are  briefly 
discussed*  as  are  applications  to  inventory  manasement* 
estimation  of  time-varyinS  mean  functions*  and  multiple 
regression. 
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I.  INTRODUCTION  AND  SUMMARY 


A.  DACKGRQUND 

The  Kalman  filter  is  3  recursive  Bayesian  least-sauares 
estimator  of  an  n-dimensional  system  state  vector  based  on  an 
m-dimensional  measurement  vector*  The  filter  may  operate  in 
a  J-dimensional  coordinate  system  where  J£mi>  J:£n*  The  basic 
assumption  is  that  each  dimension  of  the  coordinate  system 
varies  accor'dind  to  a  kth  order  Oauss-Markov  process*  The 
f-'^lman  Filter  was  developed  in  the  early  1960's  by  Kalman  and 
Bucy  Crefs*  1  and  23* 

/ 

The  Kalman  filter  may  be  used  to  obtain  an  optimal 
estimate  of  the  present  stater  a  prediction  of  future  states* 
and/or  smoothed  estimates  of  past  states*  The  current  state 
estimate  is  generally  used  to  determine  an  optimal  control 
input*  Future  state  estimates  are  used  to  determine  optimum 
present  policy*  Smoothed  past  state  estimates  are  used  for 
data  analysis  and  model  building*  Thus  the  potential  areas  of 
application  span  the  field  of  time  series  analysis* 

Applications  of  the*  Kalman  filter  are  numerous  and  the 
theory  is  being  continually  developed  and  extended*  An 
overview  of  the  development  of  linear  filtering  theory  and  an 
extensive  bibliography  may  be  found  in  Kailath  C ref *33*  A 
reasonably  clear  presentation  of  theory  and  applications  is 
contained  in  Gelb  C ref *43. 
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Perhaps  the  widest  and  most  successful  application  of 
Kalman  filtering  has  been  to  vehicle  tracking  and  control * 
Clark  C ref *53  has  written  a  particularly  lucid  description  of 
the  design  of  a  filter  for  an  anti-aircraft  gun  fire  control 
system  which  is  noteworthy  for  its  clarity  of  presentation  of 
the  underlying  theory*  It  is  evident  Crefs*  4  and  53  that 
the  design  process  is  heuristic*  and  requires  extensive 
testing  and  analysis  of  candidate  filter  configurations*  even 
when  the  process  is  well-understood  and  is  based  on  a  mature 
technology* 

The  Kal  .-.an  filter  has  also  been  applied*  with  varying 
degrees  of  success*  in  economic  models*  inventory  models*  and 
even  weather  models*  Considerable  difficulty  is  encountered 
in  model  building*  because  the  filter  design  requires  good 
estimates  of  the  variance  and  covariance  of  noise  sources*  as 
well  as  an  accurate  state  transition  model.  A  prior  estimate 
of  system  state  and  covariance  is  also  required*  which  is 
somewhat  less  critical  because  errors  in  the  prior  estimate 
decrease  with  time*  These  parameters  are  often  difficult  to 
determine  in  highly  random  processes  of  questionable 
stationarity* 

The  Kalman  filter  is  derived  and  designed  almost  entirely 
within  the  time  domain*  although  Clark  Cref*53  does  refer  to 
the  concept  of  filter  bandwidth*  The  Kalman  filter  is 
essentially  a  low-pass  filter  with  a  very  wide  transition 
band*  and  higher-order  filters  have  some  amplif ication  at  the 
mid  or  low-mid- freauency  range*  In  general*  the  stop  band 
does  not  completely  attenuate  high  frequencies*  This  allows 
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the  filter  to  attenuate  hidh-f reouency  noise  somewhat  while 
still  retaining  some  response  to  sudden  changes  of  state* 

B*  PURPOSE 

The  purpose  of  this  thesis  is  to  acouaint  the  reader  with 
the  Kalman  filter*  to  show  how  the  choice  of  various  filter 
parameters  affect  its  performance*  and  to  provide  design 
insight  through  analysis  in  the  freouency  domain*  The 
approach  is  tutorial*  and  the  reader  is  referred  to  some  of 
the  interesting  examples  which  may  be  found  in  the 
literature* 

C*  METHOD 

The  freouency  response  of  several  simple  filter  designs 
were  investigated  using  the  Fast  Fourier  Transform  program  in 
the  APL  Library  2*  The  computer  results  were  Justified 
analytically  for  the  simplest  design*  a  scalar  single-state 
filter*  Derivations  are  presented  in  appendix  A* 

D*  LEVEL  OF  PRESENTATION 

Full  understanding  of  the  theory  reouires  a  knowledge  of 
stochastic  processes  that  evolve  over  time*  as  well  as  an 
understanding  of  digital  signal  theory  in  the  freouency 
domain*  The  Fourier  transform  is  a  basic  tool*  A  full 
exposition  of  the  underlying  theory  is  clearly  beyond  the 


scope  of  this  presentation*  The  reader  is  directed  to  Larson 
and  Shubert  C ref *63  for  the  theory  of  stochastic  processes 
and  to  Hamming  Cref*73  for  the  theory  of  digital  filtering. 
As  previously  mentioned)'  Gelh  Cref*43  and  Clark  Cref*53  are 
good  references  for  the  Kalman  filter*  Bloomfield  Cref*83 
and  Brillinger  Cref.93  are  also  applicable  references*  Brown 
Cref*103  and  Box  and  Jenkins  Cref*113  contain  related 
material* 

There  ar.'  few  readers  who  are  entirely  conversant  with 
both  the  freeuency  domain  and  time  domain  approach  to  time 
series  analysis*  Nevertheless*  a  duality  exists  between  the 
two*  and  a  summary  of  the  theory  is  presented* 

Illustrative  examples  will  often  be  based  on  tracking 
models*  because  this  is  presently  the  widest  area  of 
application  of  Kalman  filters*  and  because  most  readers  will 
find  the  concepts  of  position*  velocity*  and  acceleration 
easy  to  understand*  The  concepts  are  easily  extendable  to 
other  areas*  For  example*  the  economist  may  wish  to  replace 
•velocity"  with  •trend* ♦ 

E.  SUMMARY  OF  RESULTS 

The  steady-state  gain*  bandwidth*  and  sensitivity  of  the 
linear  discrete  Kalman  filter  are  shown  to  be  completely 
determined  by  the  choice  of  the  process  and  measurement  noise 
covariances.  Filter  performance  on  stationary  or 
nearly-stationary  data  can  be  predicted  by  comparing  the 
freauency  response  of  the  proposed  filter  with  a  spectral 
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analysis  of  the  data. 

The 

wide 

transition  band 

of 

the 

amplitude  response  of 

the 

scalar 

Kalman  filter 

can 

be 

sharpened  by  multiple  passes  of  the  data  through  a 
hidher-dain  filter.  This  can  be  accomplished  simply  and 
recursively.  The  superiority  of  symmetric  smoothing  filters 
over  non-symmetric  filters  was  demonstrated.  When  used  as  a 
smoother  (by  usind  both  forward  and  backward  passes)  the 
Kalman  filter  was  as  effective  as  a  non-recursive  Gaussian 
filter.  Hidher-order  filters  were  shown  to  have  hidher 
bandwidth  and  amplif ication  as  the  order  of  the  filter  was 
increased.  A  freouency  domain  approach  to  filter  desidn  may 
provide  additional  insidht  and  enable  the  designer  to  achieve 
better  filter  performance*  particularly  when  the  system  state 
transition  model  and  noise  covariance  models  are  not  well 


understood. 


II 


THEORY 


A.  STOCHASTIC  PROCESSES  AND  STATIQNARITY 


A  continuous  stochastic  process  X<t)  is  a  Gaussian 
Process  if  t.ne  probability  densities  of  all  orders  are 
multivariate  Gaussian  densities*  It  is  a  kth  order 
Gauss-Markov  process  if  the  state  at  time  t  depends  only  on  k 
earlier  states*  If  we  should  expand  the  state  space  to  k 
states*  which  include  all  derivatives  up  to  the  (k-l)th*  the 
future  system  state  vector  will  depend  only  on  the  present 
state*  For  example*  if  the  acceleration  of  a  vehicle  is  a 
first-order  Gauss-Markov  process*  then  the  position  of  the 
vehicle  is  a  third-order  Gauss-Markov  process*  However*  if 
our  state  space  includes  acceleration  and  velocity  as  well  as 
position*  the  future  state  of  the  system  is  independent  of 
all  but  the  present  state*  If  the  random  acceleration  has 
zero  mean*  and  variance  one  over  one  time  increment*  the 
acceleration  is  a  standard  Wiener  process  W(t>*  The 
derivative  of  the  Wiener  process*  written  dW(t>*  has  zero 
mean*  unit  variance*  and  is  called  white  Gaussian  noise* 
which  may  be  thought  of  as  a  *zero-th  order*  Gauss-Markov 
process  C ref *43* 

The  standard  Wiener  process  is  not  stationary*  because 
the  variance  grows  linearly  with  time*  That  is*  the  estimate 
of  a  future  state  based  on  the  present  state  has  variance 
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that  is  a  linear  function  of 


time*  However*  the  standard 


Wiener  process  has  stationary*  independent  increments.  That 
is*  the  variance  at  time  (t+1)  given  the  state  at  time  (t)  is 
constant  and  independent  of  t. 

A  stochastic  process  X<t>  is  wide-sense  stationary  if  and 
only  if  it  has  a  constant  mean  function*  and  a  correlation 
function  such  that  C ref *61 

R^t.t  s*tz+  s>  =  R  <t,  *t,>  =  R  <tx-t,> 

that  is*  the  correlation  function  of  the  process  is 
independent  of  an  arbitrary  time  shift  s*  A  Gaussian  process 
is  strictly  stationary  if  and  only  if  it  is  wide-sense 
stationary  C ref. 63* 

The  Gauss-Markov  assumption  makes  possible  the 
development  of  theory  and  applications*  because*  in  general* 
any  linear  operation  performed  on  a  Gaussian  process  results 
in  another  Gaussian  process*  and  the  Markov  property  allows 
consideration  of  only  the  present  state*  disregarding  all 
previous  states* 

B.  THE  PHILOSOPHICAL  CONCEPT  OF  STATIONARITY 


A  f reouency-domain  analysis  of  a  stochastic  process  is 
only  meaningful  if  the  process  is  stationary*  If  the  process 
were  changing  over  time*  the  spectrum  would  change  over  time* 
Since  the  spectrum  can  only  be  analyzed  by  means  of  data 
taken  over  time*  such  analysis  of  a  nan-stationary  process 
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would  be  meaningless*  However;  if  the  process  is 
•ouasi-stationary • ;  that  is;  it  exhibits  stationary 
statistics  for  a  while ;  then  undergoes  a  change;  then  settles 
down  to  stationarity  again;  the  freeuency  approach  is  still 
useful ;  although  inaccurate  over  the  transition  period*  As 
an  example ;  consider  an  airplane  subject  to  random 
accelerations  due  to  air  turbulence*  An  appropriate  model 
might  be  a  third-order  Gauss-Markov  process  as  long  as  the 
airplane  maintains  a  straight  path  or  turns  at  a  constant 
acceleration*  However;  the  pilot's  inputs  to  initiate  or 
terminate  a  maneuver  would  result  in  brief  periods  of 
non-stationarity*  and  the  model  would  perform  inadeouately 
during  and  immediately  after  the  transition  period* 

It  may  be  argued  that  every  practical  process  can  be 
considered  stationary  over  infinite  time.  If  the  process  is 
random;  it  represents  an  ensemble  of  possible  paths;  of  which 
any  realization  in  terms  of  real-World  data  is  only  one 
possible  path;  and  may  or  may  not  be  closely  representative 
of  the  ensemble*  Uhen  dealing  with  reality;  we  are  often 
forced  to  assume  stationarity  in  order  to  make  analysis 
possible;  and  often  we  obtain  good  results  even  though  we  can 
never  know  whether  or  not  the  assumption  of  stationarity  is 
really  valid* 

C.  DUALITY  OF  THE  TIME  AND  FREQUENCY  DOMAINS 


1 ♦  Fourier  Series 

A  very  wide  class  of  mathematical  functions  may  be 
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represented  by  the  Fourier  series  Cref.122  as  follows! 


s<t>  =  a#  +  ^(^cos  nt  +  bfl  sin  nt) 


where  0  £  t  :£  2lf* 

Existence  and  convergence  of  this  series  representation 
reauire  only  that  g<t)  be  everywhere  single-valued*  and 
possess  a  finite  number  of  maximal*  minimal*  and  finite 
discontinuities*  The  function  g<t>  need  not  be 
differentiable*  Any  function  meeting  the  above  criteria  can 
be  thought  of  as  a  constant  mean  function  ae*  plus  an 
infinite  series  of  sines  and  cosines  of  integral  freouencies 
and  various  amplitudes.  Of  courser  the  independent  variable 
t  must  be  shifted  and  scaled  to  the  interval  C  0r27TT*  Note 
that  the  lowest  freouency  present*  aside  from  the 
zero-f reouency  mean*  is  one  cycle  for  the  span  of  g(t). 
Among  the  functions  meeting  the  criteria  are  a  seuare  pulse* 
an  impulse*  and  any  manifestation  of  a  random  walk*  In 
practice*  the  Fourier  analysis  of  a  function  s(t)  reouires 
the  truncation  of  the  infinite  Fourier  series*  This  results 
in  a  smooth  least-souares  approximation  to  the  function  g(t>* 
There  are  ripples  in  the  approximation  if  the  function  s(t) 
is  not  differentiable  or  if  the  truncation  is  too  severe* 
This  is  known  as  the  Gibbs  phenomenon*  and  is  illustrated  in 
figure  1*  which  was  taken  from  Hamming  C ref *73*  By  taking  a 
sufficient  number  of  terms  in  the  Fourier  expansion*  we  can 
improve  the  closeness  of  the  approximation* 
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2*  Basic  Concept 

The  basic  concept  of  the  duality  between  the  time  and 
freauency  domains  is  so  simple  that  it  often  gets  lost  in  3 
forest  of  Fourier  transforms*  The  time  period  is  the 
reciprocal  of  frequency*  The  basic  relationship  is 

v/27T  *  f  -  1/T 

where  v  is  the  freauency  in  radians/unit  timer  f  is  the 
freauency  in  cycles/unit  timer  and  T  is  the  time  period  for 
one  cycle.  Stated  simply r  freauency  is  the  inverse  of  the 
time  period. 

3.  Discrete  Data  and  the  Sampling  Theorem 

The  digital  computer  allows  the  efficient  analysis  of 
continuous  phenomena  by  means  of  discrete  approximations.  We 
saw  earlier  that  the  lowest  freauency  contained  in  a  Fourier 
expansion  of  a  function  g<t>  was  the  reciprocal  of  the  time 
span  covered  by  the  function.  Similarlyr  the  famous  Sampling 
Theorem  Crefs.  6  and  73  states  that  if  a  function  S<t>  in 
continuous  time  is  sampled  at  constantr  discrete  time 
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intervals  kt  (that  is*  at  a  rate  of  i/kt>*  then  the  highest 
observable  freouencs  is  0*5  cades  per  measurement  interval 
kt*  This  means  that  at  least  two  observations  are  reouired  in 
each  cacle  in  order  to  observe  that  particular  freouenca* 
The  freouenca  0.5/kt*  usualla  written  simpla  0*5*  is 
referred  to  as  the  Naouist  freouenca.  The  result  is  the 
aliasing  phenomenon*  which  is  familiar  to  most  moviegoers* 
During  the  chaser  the  stagecoach  wheels  appear  to  stop  or 
rotate  slowla  backwards  when  the  rate  of  rotation  of  the 
wheel  spokes  (spokes/sec)  exceeds  1/2  the  camera  rate 
(f rames/sec) ♦  When  higher  freouencies  exist  in  the  function 
g(t)  sampled  at  a  rate  kt*  thea  are  folded  back  and  appear  in 
the  freouenca  spectrum  of  the  sampled  data  as  freouencies 
less  than  the  Naouist  freouenca*  The  sampling  theorem  shows 
that  a  spectral  analasis  of  discrete  data  is  onla  meaningful 
over  the  Naouist  interval  C  -0*5/bt*  +0.5/l^t3« 

4*  The  Discrete  Fourier  Transform 

Ana  function  g(t>  for  which  a  convergent  Fourier 
series  exists  mas  be  represented  in  the  freouenca  domain  in 
terms  of  real  and  imadinars  parts*  or  in  terms  of  amplitude 
and  phase  angle*  as  a  function  of  freouenca*  It  should  be 
noted  that  the  function  d(t)  mas  also  be  complex-valued*  but 
we  will  deal  with  onla  real-valued  functions*  In  the 
continuous  domain*  the  formulas 

g(t)  » 

and 


si' 


G(v)  exp(ivt)  dv 
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G(v) 


exp(-ivt)  dt 


represent  3  Fourier  transform  pair*  The  freeuency  response 
G<v>  completely  determines  the  time  function  s(t)  and 
conversely* 

If  the  function  si<t)  is  sampled  at  intervals  t  « 
0»  1  r2f  ♦ » •  »n  the  time-to-f  reetuency  transformation  becomes 


A 

G(v)  -  21d<t)  exp(-ivt>  t  =  0»l»*.*»n 

t«o 

which  is  defined  only  on  the  Nyouist  interval  C-7r,7T;],  here 
defined  in  radians*  The  formula  may  be  written  in  a  more 
'familiar  form  by  usind  the  Euler  relation 


exp<-ivt)  =  cos  vt  -  i  sin  vt 


as 


It 

G<v)  =  2Jd<t)  (cos  vt  -  i  sin  vt) 

t*o 

which  is  continuous  in  v  on  the  interval  C-irrTTl*  The  Fourier 
transform  is  a  bit  difficult  to  handle  analytically  for  all 
but  the  simplest  functions t  but  the  discrete  Fourier 
transform  is  Generally  easy  to  compute  by  use  of  a  Fast 
Fourier  Transform  <FFT>  prodram  available  in  most  computer 
libraries*  The  output  will  generally  be  a  very  close  discrete 
approximation  to  G<v)r  if  the  span  of  d<t>  is  larste  enoudh* 
The  inverse  transformation  can  also  be  made* 
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Since  G<v>  is  complex  valued  whenever  the  function  S<t> 
is  not  symmetric*  it  is  often  useful  to  represent  it  in  terms 
of  amplitude  and  phase*  The  amplitude  is 


|G(v> 


G(v)G(-v) 


-I 


CRe<vn*+  CIm<v)3l 


where  Re(v)  and  Iro(v)  are  the  real  and  imaginary  parts  of 
G(v)*  The  phase  angle  is 


8<v)  =  arctan  Im(v)  /  Re(v) 


D.  THE  DOOB-MEYER-FISK  DECOMPOSITION 


In  most  practical  applications*  a  finite-variance 
sample-continuous  stochastic  process  X(t>  can  be  written 


T  r 

X(T)  =  XCO)  +  /a<  t)dt  +  yB(t)dW(t) 

o 


where  X(0)  is  the  initial  value  of  that  process*  A(t)dt  is 
predictable*  smooth  behavior  determined  by  a  set  of 
deterministic  differential  eauations  describing  the  system* 
and  B(t)dW(t)  is  noise*  where  dU(t>  is  white  Gaussian  noise* 
and  B<t)  is  a  smooth  transformation  that  is  sometimes  thought 
of  as  "coloring*  the  noise*  Such  a  representation  is  called 
the  Doob-Meyer-Fisk  decomposition  C ref. 63*  which  may  be 
thought  of  as  separating  the  process  into  a  signal  and  noise* 


Several  important  points  must  be  made  with  regard  to  this 


eouation.  It  is  not  intended  here  that  the  expression  be 
evaluated  analytically*  The  integral  B<t)dW<t)  is  an  IttJ 
intedral*  which  is  not  even  a  stochastic  version  of  a 
StieltJes  intedral  Cref.Al.  Also*  although  the  processes 
A<t)  and  B<t)  are  smooth  functions  that  may  be  considered 
deterministic  representations  of  system  behavior*  they  are 
not  necessarily  known  to  the  observer*  even  when  an  adeeuate 
technological  representation  exists. 

Consider  adain  our  piloted  aircraft  beind  tracked  by  a 
radar.  The  process  A<t)  represents  the  dynamics  of  the 
airframe*  as  affected  by  the  control  inputs  of  the  pilot* 
which  are  unknown  to  the  radar  observer.  The  process  B(t) 
consists  of  several  parts.  One  is  the  measurement  process* 
which  may  or  may  not  be  known  to  the  radar  observer.  For 
example*  a  nutatinst  radar  antenna  midht  impose  some  periodic 
error  in  the  measurement*  which  would  be  manifested  in  the 
process  B<t).  Overlaid  on  this  midht  be  a  white  Gaussian 
noise  measurement  error.  Air  turbulence  could  also  be 
represented  as  white  Gaussian  noise*  which*  however*  could 
only  be  manifested  throudh  deterministic  airplane  dynamics. 
There  are  those  who  would  ardue  that  the  pilot  should  also  be 
modelled  as  a  random  variable.  In  any  event*  the  process 
B(t)  midht  be  further  decomposed  into  several  processes*  here 
at  least  airframe  response  to  air  turbulence  and  periodic 
radar  antenna  dynamics. 

The  vital  observation  is  that  if  the  frequency  content  of 
the  processes  A<t)  and  B(t)  are  known  to  be  different*  they 
can  be  partially  separated  by  a  spectral  analysis  of  the 


data 


In  our  example*  aircraft  have  natural  dynamic  response 


frequencies  in  all  control  axes.  These  can  bo  estimated 
closely*  even  for  enemy  airplanes*  and  are  Generally  similar 
amend  similar  types  of  airplanes*  although  they  vary  with 
airspeed.  It  is  physically  impossible  for  the  airplane  to 
respond  faster  than  its  highest  natural  dynamic  frequencies.* 
Any  frequency  content  higher  than  this  must  be  noise.  If  the 
radar  system  dynamics  are  of  a  higher  frequency  than  this* 
they  can  also  be  separated.  The  pilot  will  take  advantage  of 
the  full  response  rate  of  the  airplane  only  very  rarely. 
Therefore*  low  freouency  components  are  most  likely  due  to 

pilot  maneuvers.  Of  course*  since  white  Gaussian  noise  has 
a  flat  freouency  spectrum  as  a  result  of  aliasing  Cref.73*  it 
is  impossible  to  separate  all  of  the  noise  from  the  signal. 
However*  it  is  often  possible  to  remove  quite  ?.  bit  of  it. 

E.  DIGITAL  FILTERS 

A  digital  filter  is  a  linear  transformation  applied 
iteratively  to  a  set  of  data  points.  The  purpose  here  is  to 
separate  noise  from  the  signal.  The  simplest  digital  filter 
is  the  simple  average*  which  estimates  the  mean  value  from 
the  data*  and  smooths  out  all  fluctuations.  The  most  general 
form  of  the  digital  filter  was  stated  by  Hamming  Cref.71  as 

•• 

x ( t )  =  ^a(k)  z(t-k)  +  b ( k )  x(t-k) 

4»J 

where  the  estimate  x<t)  at  some  point  t  is  a  linear 
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combination  of  the  data 


points 


<t-k)*  and  perhaps  of  the 


previous  estimates  x(t-k>*  The  coefficients  a<k)  and  b(k) 
are  weighting  coefficients  and  may*  of  course*  be  aero*  As  a 
result  cf  the  sampling  theorem*  the  filtering  process  is 
meaningless  unless  the  measurements  z(t-k>  are  made  at 
eoually  spaced  intervals  along  the  t  axis*  where  t  is 
usually*  but  not  necessarily*  time* 


1 ♦  Some  Classif ications  of  Digital  Filters 

Digital  filters  may  be  classified  as  symmetric  or 
non-symmetric*  and  as  recursive  or  non-recursive*  A 
symmetric  non-recursive  filter  is  one  in  which  all  b(k)  eauai 
zero  and  all  a<k)  *  a<-k)*  such  as  the  filter 


:< t )  «  0*2  2<t-l>  +  0,6  z<t>  +  0.2  z(t+l>. 


An  example  of  a  recursive  filter  is 


x(t)  *  a  2<t)  +  b  x<t~l>  0<a<l *  b  =  l-a 


which  is  not  symmetric*  This  particular  filter  may  be 
expressed  as 


x(t)  =  a  2 ( t )  +  bC  a  z<t-l>  +  bC  a  2<t-2)  +  ♦♦♦♦in 


which  reduces  to 


x<t)  =  a  2<t)  +  ab  z<t-l)  +  ab*2<t-2)  +♦♦.+  abn2<t-n)  +♦♦♦ 
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The  recursive  filter  extends  to  the  infinite  past*  although 
the  coefficients  eb  will  approach  zero*  if  |bl<l*  In  this 
case*  a  recursive  filter  can  be  closely  approximated  by  * 
non-recursive  filter*  A  primary  advantage  of  the  recurrxv-“ 
filter  is  that  old  data  need  not  be  stored*  New  estimates 
may  be  computed  simply  and  rapidly  as  time  evolves*  This  is 
an  important  advantage  for  real-time  applications* 

2*  Applications  of  Digital  Filters 

Digital  filters  are  used  to  separate  a  signal  from 
noise*  to  separate  various  frequency  components  of  a  signal* 
and/or  to  perform  such  mathematical  functions  as  integration 
and  differentiation*  A  review  of  Simpson's  rule  and  the 
Trapezoidal  rule  should  convince  the  reader  that  these 
numerical  integration  techniques  are*  in  fact*  recursive 
digital  filters*  Sometimes  a  filter  has  two  purposes*  For 
example*  it  might  be  desirable*  n  estimating  velocity  from 
successive  observations  of  position*  to  simultaneously 
differentiate  and  remove  high  frequency  noise*  When  a  filter 
is  used  to  stop  part  of  the  frequency  spectrum*  it  is 
referred  to  as  a  "low-pass**  "high-pass"*  "band-pass"*  or 
"band-stop*  filter*  depending  on  its  function* 

3*  Analysis  of  Digital  Filters 

In  the  time  domain*  a  digital  filter  is  described 
completely  by  its  impulse  response  function*  which  is  nothing 
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more  than  the  response  of  the  filter  to  data  consisting  of  a 
string  of  zerosr  a  single  one*  followed  ba  zeros*  The  output 
of  the  filter  is  then  simply  the  weighting  coefficients  a(k>* 
If  the  filter  is  recursive*  we  might  not  be  able  to  deduce 
the  recursive  form  from  the  coefficients  a<k)»  but  that  will 
not  concern  us  here*  The  Fourier  transform  of  the 
impulse-response  function 

N 

H<v)  =  5”a<k)  exp(-ivk) 

will  completely  specify  the  frequency  response  of  the  filter* 
If  the  filter  is  symmetric*  there  will  be  no  imaginary  part* 
and  hence  no  phase  shift*  If  the  filter  is  recursive*  it 
cannot  practically  be  symmetric*  and  the  summation  will 
generally  run  from  zero  to  infinity.  That  is*  the  impulse 
response  will  extend  infinitely  far  into  the  future*  which 
means  that  the  filter  remembers  all  of  the  past* 

The  duality  of  the  time  and  freauency  domains  allows  us 
to  specify  a  desired  freauency  response  and  to  design  an 
appropriate  filter  by  calculating  filter  weights*  or  to 
analyze  an  existing  filter  by  calculating  the  freauency 
response  from  the  filter  weighting  coefficients* 

F*  DATA  ANALYSIS  AND  EXPERIMENT  DESIGN 

No  digital  filter  should  be  applied  to  data  analysis 
without  a  clear  idea  of  the  effect  of  the  filter  upon  the 
data*  Slutsky  and  Yule  first  noted  that  some  smoothing 
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formulas  induced  periodic  functions  in  the  smoothed  estimate 
that  were  more  the  effect  of  smoothing  than  of  the  original 
data  C ref. 73,  A  spectral  analysis  of  representative  raw  data 
can  be  helpful  in  deciding  on  an  appropriate  filtering 
technique*  However *  such  data  as  economic  time  series  or 
weather  data  are  typically  very  noisy*  are  based  on  a 
relatively  short  run  of  data*  and  cannot  be  described  by  an 
adeouate  technological  model*  The  analyst  must  be  aware  of 
these  problems.  Sometimes  there  are  no  good  solutions*  but  a 
spectral  representation  may  produce  freouencies  that  can  be 
explained  on  rational  grounds* 

Another  potential  pitfall  is  a  result  of  the  sampling 
theorem*  Consider  the  timely  example  of  an  air  pollution 
model.  It  would  be  reaonable  to  suspect  that  air  pollution 
would  follow  at  least  a  daily  cycle*  or  perhaps  an  eight  hour 
cycle  if  morning  and  evening  rush  hours  were  considered* 
Daily  samples  of  air  pollution  could  not  hope  to  uncover 
cycles  of  a  shorter  period  than  ?very  two  days*  Samples 
every  four  hours  would  be  marginally  adeeuate*  Hourly 
samples  would  be  necessary  for  a  good  analysis*  Additionally* 
recall  the  reouirement  for  eaually-spaced  sampling  intervals* 
For  various  reasons*  the  analyst  may  have  no  control  over 
data  collection*  However*  he  must  always  understand  what  has 
been  done*  or  could  have  been  done*  to  the  data*  as  well  as 
what  he  is  doing  to  it*  in  order  to  avoid  erroneous 
conclusions* 
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III.  THE  LINEAR  DISCRETE  KALMAN  FILTER 


A.  DESCRIPTION 

The  linear  discrete  Kalman  filter  is  a  recursive  Bayesian 
least-sauares  estimator  of  the  state  vector  of  a  linear 
system  based  on  a  vector  of  noisy  measurements  made  at 
discrete  time  intervals*  The  process  to  be  estimated  is 
assumed  to  be  an  n-state  Gauss-Markov  process  of  order  k t 
subject  to  process  noise  W  with  aero  mean  and  covariance 
matrix  Q.  The  process  is  observed  by  an  m-dimensional 
measurement,  subject  to  measurement  noise  V  (not  to  be 
confused  with  freouency  <v)>  with  zero  mean  and  covariance 
matrix  R.  The  filter  reouires  a  prior  Bayesian  estimate  of 
system  state  and  covariance*  The  recursive  estimate  of  system 
state  at  time  t  is  obtained  by  the  formula 

X(t!t)=  X< t ! t-1 )  +  K(t ) CZ( t )-HX( t I t-1  )3 

where 

X(t't)  state  estimate  based  on  current  measurement 

Z(t)  current  measurement 

X(t't-l)  state  estimate  prior  to  current  measurement 

K(t)  Kalman  slain  matrix  (to  be  discussed  later) 

H  observation  matrix,  which  is  constant 

The  derivation  of  the  Kalman  filter  equations  may  be 


found  in  Gelb  C ref *43*  A  summary  of  the  filter  equations  is 
presented  in  figure  2 *  which  should  be  consulted  in  order  to 
follow  the  subsequent  discussion* 

In  general *  the  state  model  represents  a  dynamic  system* 
that  is*  one  which  changes  with  time*  The  extrapolation  of 
the  state  estimate  to  the  time  of  tb’  next  observation  is 
obtained  by  the  formula* 

X<t+1 it)  =  $  X<tlt) 

where  $  is  the  state  transition  matrix.  The  observation 
process  occurs  according  to  the  conceptual  relation 

Z<t>  *  HX(t)  +  V 

where  X(t)  is  the  true  system  state,  observed  through  the 
observation  matrix  H*  and  V  represents  measurement  noise* 
which  is  assumed  to  be  a  Gaussian  random  variable  with  zero 
mean  and  covariance  matrix  R*  Note  that  the  process 
represented  by  this  formula  is  assumed  to  occur  in  the  real 
world*  The  computation  does  not  occur  in  the  filter*  Rather* 
the  measurement  Z(t)  is  an  input  to  the  filter* 

In  the  linear  Kalman  filter*  the  gain  K(t>  does  not 
depend  in  any  way  on  the  data*  It  depends  only  on  the  model* 
and  is  therefore  extremely  sensitive  to  assumptions*  Gain  is 
calculated  according  to  the  formula 

T  t  -1 

K(t)  =  P<t)H  CHP<t>H  +  R1 
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where  P(t)  is  the  covariance  in  the  system  state  estimate 
prior  to  the  current  measurement  and  R  is  the  covariance  of 
the  measurement  error*  The  covariance  is  updated  according 
to  the  formula 

£<t>  =  CI-K<t)H3  P<t> 

where  £<t>  is  the  state  covariance  given  the  current 
measurement)!  and  I  is  the  identity  matrix*  The  covariance  is 
extrapolated  to  the  time  immediately  prior  to  the  next 
observation  by  the  formula 

P<t+1>  =  J  Z(t)  JT+  Q 

where  $  is  the  state  transition  matrix  and  Q  is  the 
covariance  of  the  process  noise*  Combining  the  above  two 
eouations  shows  that  the  covariance  of  the  state  estimate  at 
the  time  of  the  current  measurement  depends  on  the  previous 
covariance  according  to  the  formula 

£<t>  =  CI-K<t)H3  C$  Z<t-1>  $T +  Q3 

Filter  performance  is  very  dependent  on  adeouate 
modelling*  particularly  on  the  state  transition  model  $  and 
the  choice  of  noise  covariances  R  and  Q*  To  a  lesser  extent* 
performance  also  depends  on  the  initial  estimates  of  system 
state  X<1!0>  and  covariance  P<1!0>*  However*  the  latter 
parameters  are  less  important  because  their  effects  decrease 
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with  time*  If  the  matrices  R>  G*  5*  and  H  are  constant  in 
time*  the  Sain  K(t)  and  covariance  matrices  £(t)  and  F’(t) 
eventually  reach  a  steady  state*  and  are  completely 
determined  by  R*  G*  jg*  and  H* 

For  a  given  linear  filter*  it  will  be  shown  that  filter 
gain*  covariance*  and  freouency  response  will  be  completely 
determined  by  the  choice  of  R  and  G. 

B.  THE  SCALAR  KALMAN  FILTER 

The  multi-state  Kalman  filter  is  a  powerful  computational 
device*  However*  it  is  difficult  and  often  impossible  to 
manipulate  in  closed  form  because  of  the  freauent  occurence 
of  singular  matrices*  An  analysis  of  a  single-state  (scalar) 
filter  can  be  used  to  illustrate  the  mechanics  of  the  Kalman 
filter*  and  to  aid  in  developing  an  intu'tive  understanding* 
In  the  discussion  that  follows*  it  is  assumed  that  all 
matrices  are  scalars*  and*  in  particular*  Q  and  H  eeual  one* 
Matrix  notation  is  preserved  for  clarity*  Derivations  may  be 
found  in  appendix  A* 

1  *  Transient  and  Steady-State  Gain 

It  can  be  shown  (appendix  A)  that  the  scalar  Kalman 
gain  can  be  expressed  recursively  as 

K(t )  =  K(t-l)  +  Q/R 

K(t-l)  +  G/R  +  1 

When  the  filter  reaches  steady-state*  the  gain  is  constant 


3nd 


K  =  +  /j3*  +  Q 

2R  J 4R*  R 

The  inverse  relationship  is 

Q  =  j£ 

R  1-K 

Thus*  the  variance  ratio  Q/R»  which  is  the  ratio  of  process 
noise  variance  to  measurement  noise  variancef  completely 
determines  the  steady-state  Sain*  The  steady-state  filter  is 
completely  described  by  the  formula 

X<t>  =  K  Z<t>  +  < 1-K)  X(t-l) 


2»  Frequency  Response 

Lettind  K  =  a  and  <1-K>  =  b»  the  impulse-response 

function  G(t)  may  be  written 


G ( t )  =  ab  f  t=0 r 1 r  2  f 
The  Fourier  Transform  is 


-J 


H(v)  =  I  G(t)  e;;p  <-ivt)dt 


H(v)  »  a  ^Cb  eMp(-iv)]^ 

t«o 


H(v)  =  a  /  Cl  -  b  e*p<-iv)3 


Since  the  filter  is  not  symmetric ?  the  freauency  response 
H(v)  has  both  reel  end  imaginary  parts*  The  amplitude  may  be 
written 

A  =  |H(v)|  =  yH<v)H<-v)  =  a  /  J 1  +  b*  -2b  cos  v 

which  reduces  to  ^ 

A  =  /  Q/R 

J  Q/R  +  2(1  -cos  v) 

The  phase  angle  may  be  written 

8(v)  =  arctan  /  -b  sin  v  j 

ll  ~b  cos  v/ 

The  angle  for  maximum  phase  shift  is 

v(max  8)  «  arccos  b  =  arccos  <1-K) 

v  (max  8)  =  arccos  <1  +  Q/2R  -  Jat/4RZ  + Q/R  ) 

Therefore?  the  variance  ratio  Q/R  also  completely  specifies 
the  steady-state  freauency  response  of  the  filter*  Amplitude 
and  phase  relationship  for  several  values  of  gain  are  plotted 
in  figures  3  and  4* 

It  is  evident  that  high  Q/R  (high  gain)  reduces  the  phase 
lag  of  the  filter  but  allows  more  of  the  hish-f reouency 
components  to  pass*  Conversely?  low  Q/R  (low  gain) 
attenuates  more  of  the  high-f reouency  components?  at  the 
expense  of  an  increased  phase  lag*  Notp  that  even  at  very 
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FREQUENCY  RESPONSE  OF  SCALAR  KALMAN  FILTER 
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low  gain  <low  G/R>  not  all  of  the  hish-f reouency  component  is 


attenuated*  and  the  phase  lag  is  ouite  severe*  The  slope  of 
the  amplitude  change  is  ouite  shallow*  implying  that 
attenuation  increases  gradually  as  freouency  increases*  This 
is  a  conseouence  of  the  assumptions  and  performance  will  not 
he  adeouate  if  the  data  does  not  represent  a  Gauss-Markov 
process*  but  in  fact  represents  some  phenomena  changing  with 
time* 

C.  IMPROVING  THE  FILTER 

The  transition  band  of  the  filter  can  be  sharpened*  and 
more  of  the  high-  freouency  components  eliminated*  by  running 
the  data  through  two  filters  in  series*  The  basic  scalar 
filter  was 


x<t)  =  a  z<t>  +  b  x(t-l ) 

where  a  »  K  and  b  =  <1-K>.  Running  the  data  through  the 
filter  again*  we  obtain  a  new  estimate  y<t>  *  where 

y(t>  =  a  x<t>  +  b  y<t-i) 

It  should  be  evident  that  we  can  accomplish  this  3ll  in  one 
step  as 

y<t>  =  a*s<t)  +  2b  y(t-l)  -  b*  y<t-2) 


We  need  only  to  save  one  additional  previous  estimate  y(t-2) 
as  well  as  y(t~l>.  The  impulse  response  function  is 


g(t )  =  (t+l)a*  fa*  »  t  =  0*1*... 


We  have  performed  a  convolution  in  the  time  domain*  which 
corresponds  to  a  multiplication  in  the  freouency  domain. 
This  may  not  be  exactly  what  we  want.  Let  us  suppose  that  we 
want  the  weighting  coefficient  for  the  present  data  point 
z(t)  to  be  0.27  in  both  cases.  This  reouires  a  =  0.27  for 
the  basic  filter  and  a  =  ^0.27  =  0.52  for  the  double  filter. 
The  impulse  response  function  for  both  filters  is  presented 
in  figure  5.  Note  that  the  double  filter  forgets  the  past 
more  readily.  The  amplitude  and  phase  shift  for  both  filters 
is  presented  in  figures  6  and  7.  The  sain  for  the  scalar 

I 

filter  was  0.27*  corresponding  to  a  variance  ratio  <Q/R)  of 
0.1.  Also  (note  that  the  double  filter  has  somewhat  better 
hish-f reouency  attenuation*  somewhat  less  attenuation  at  low 
freouencies*  and  a  slightly  sharper  (steeper)  transition.  We 
would  therefore  expect  it  to  be  a  bit  better  at  separating  a 
low-f reouency  signal  from  noise.  At  low  freouencies*  the 
double  filter  has  less  phase  shift.  However*  phase  shift  is 
more  severe  at  freouencies  above  0.5  radians. 


D.  PERFORMANCE  COMPARISON 


The  basic  scalar  filter  and  the  improved  (double)  filter 


37 
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were  compared  using  the  data  illustrated  in  figure  8*  A 
series  of  independent  normal  random  numbers  with  zero  mean 
and  variance  ten  were  generated*  to  simulate  measurement 
noise*  They  were  added  to  a  nominal  function  which  was  a 
combination  of  steps  and  a  ramp*  The  ramp  function  rises  1 
unit  each,  measurement  interval*  which  corresponds  to  one 
standard  deviation  in  process  noise  G*  so  at  least  during  the 
ramp  function*  the  data  corresponds  to  the  filter  design 
variance  ratio  Q/R  =  0*1.  A  spectral  analysis  of  the  data  is 
presented  in  figure  9.  The  signal  and  noise  are  presented 
separately  and  in  combination*  Note  that  the  noise-only  scale 
is  expanded*  The  noise  spectrum  is  irregular*  but  overall 
Quite  flat.  The  signal  consists  mostly  of  very  low 
freouencies*  but  also  has  some  high  freauencies*  This  would 
be  expected*  since  step  and  ramp  functions  reauire  very  high 
freauencies  in  their  Fourier  expansion*  The  hish-f reauency 
signal  is  submerged  in  noise.  The  data*  of  course*  does  not 
fulfil  the  assumptions  from  which  the  Kalman  filter  is 
derived.  However*  the  real  world  seldom  does  either*  We  are 
looking  for  robustness* 

The  filters  were  first  tested  on  the  signal  alone*  The 
results  are  presented  in  figure  10*  It  can  be  seen  that 
neither  filter  can  respond  instantaneously  to  the 
discontinuities  in  the  function*  since  high  freauencies  are 
attenuated*  Both  filters  lag  after  discontinuities  and 
during  the  ramp  rise*  This  is  a  conseauence  of  the  non- 
symmetric  nature  of  the  filters  and  illustrates  the  phase 
lag*  Note  that  the  double  filter  performs  a  bit  better* 
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remaining  closer  to  the  data  throughout  (it  is  a  peculiarity 

of  the  plotting  routine  that  only  one  symbol  will  be  plotted 

when  ever  data  points  coincide*  Therefore  *  data  points  that 

do  not  appear  should  be  regarded  as  occurring  simultaneously 

with  the  ones  that  do  appear)* 

The  filters  were  then  tested  on  the  noisy  data*  The 

results  are  presented  in  figures  11  and  12*  The  latter  plot 

has  the  data  suppressed  so  that  the  scale  can  be  increased 

and  more  resolution  obtained* 

The  same  trends  can  be  observed  as  were  previously*  The 

double  filter  lags  less  during  the  trend  and  transitions* 

1 

The  double  filter  appears  to  follow  the  noise  a  bit  more 

closelyr  but  overall  it  follows  the  signal  better  than  the 

basic  filter*  The  average  variance  between  the  signal  and  the 

filter  was  4.20  for  the  double  filter  and  5*10  for  the  basic 

filter*  which  was  an  18%  improvement  for  this  simple 

1 

modification*  The  improvement  is  due  to  the  fact  that  the 
double  filter  weights  more  recent  ,  data  more  heavily*  and 
remembers  less  of  the  past  than  the  scalar  Kalman  filter* 
even  though  the  weight  on  the  present  observation  is  the 
same* 

This  simple  experiment  is  only  intended  to  acquaint  the 
reader  with  possible  improvements  to  the  Kalman  filter*  Like 
any  tool*  the  Kalman  filter  should  not  be  applied 
indiscriminantly ♦  The  interested  reader  is  referred  to 
Hamming  Cref*73  for  the  basic  principles  of  digital  filter 
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FIGURE  11 


E.  THE  TRANSIENT  CASE 


The  sophisticated  reader  has  no  doubt  noticed  that  the 
steady-state  scalar  Kalman  filter  is  eauivalent  to  the 
Box-Jenkins  IMA<0»1>1)  model  Cref.113  and  to  Brown's 
exponential  smoothing  model  Cref.103.  Zehna  criticised  the 
exponential  smoothing  model  Href. 133*  noting  the  bias  would 
occur  if  the  steady-state  model  was  applied  with  an 
inappropriate  prior  estimate.  Bessler  and  Zehna  C ref. 143 
developed  a  gain  schedule  which  they  call  finite  exponential 
smoothing.  Their  formula  for  gain  is 

a<t)  *  a/ (l-b*) 

where  a  is  the  steady  state  gain*  b  =  (l-a>»  and  a<t>  is  the 

« 

gain  schedule  as  a  function  of  time.  It  is  similar  to  the 
Kalman  gain  schedule  if  the  initial  Kalman  gain  K<0>  is  t 
chosen  as  one.  In  both  models*  no  prior  estimate  is 
reouired.  The  weight  on  the  first  observation  is  one.  A 
comparison  of  the  two  models  is  illustrated  in  figure  13*  for 
a  steady-state  gain  of  0.2  and  an  initial  gain  of  1.  The 
Kalman  gain  was  calculated  according  to  the  recursive  formula  • 
in  section  III. A. 

The  Kalman  filter  gain  converges  faster*  although  the 
difference  is  not  great.  The  scalar  Kalman  filter  possesses 
two  other  advantages  over  the  finite  exponential  smoothing 
techniaue.  First*  if  a  good  prior  estimate  does  exist*  the. 
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FIGURE  13 


Kalman  filter  allows  the  initial  dain  to  be  chosen  as  less 
than  one*  and  its  value  will  be  determined  by  the  assumed 
covariance  of  the  prior  estimate*  Also*  the  Kalman  model 
forces  the  analyst  to  at  least  ink  about  the  concepts  of 
measurement  noise  and  process  noise*  and  to  estimate  the 
noise  variance  ratio  Q/R* 

As  noted  before*  a  freouency  analysis  of  the  transient 
case  is  not  appropriate*  However*  it  can  be  thoudht  of  as  a 
case  of  transient  bandwidth*  The  filter  is  initially  set  to  a 
hidher  dain  than  steady-state  dain*  If  there  is  no  prior 
estimate  available*  the  initial  dain  is  one*  and  the  filter 
is  initially  an  infinite-bandwidth  or  all-pass  filter*  As 
data  are  acouired  dain  drops  and  the  bandwidth  narrows  until 
steady-state  conditions  are  achieved*  The  concept  of 
transient  bandwidth  is  important  to  the  subject  of  adaptive 
filterind*  to  which  we  will  return* 

F.  HIGHER  ORDER  FILTERS 

The  main  beauty  of  the  Kalman  filter  is  not  in  its 
statistically  unbiased  method  of  calculatind  dain*  but  in  its 
powerful  matrix  formulation*  which  allows  it  to  be  applied  as 
a  multi-dimensional  model  incorporated  any  order  of 
differencind  desired*  As  the  state  space  is  increased*  it 
Quickly  becomes  impossible  to  analyse  the  filter 
analytically*  Hidh-order  multi-dimension  filters  can  also 
easily  exceed  the  capacity  of  present  didital  computers  for 
real-time  applications*  Fortunately*  it  has  been  found  that 
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the  state  space  can  be  reduced  and  dimensions  decoupled  with 


very  little  degradation  in  the  overall  accuracy  of  the  state 
estimate  Crefs.  4  and  53*  For  example*  if  a  12-state  model 
can  be  reduced  to  9  states  and  can  be  adeouately  represented 
by  three  3-state  models*  the  matrix  calculations  can  be 
considerably  simplified  and  speeded  up. 

We  will  examine  a  second-order  (first  difference) 
filter*  which  can  be  used  to  estimate  trend*  or  velocity.  We 
will  use  the  latter  term.  Position  and  velocity  are  to  be 
estimated  bosed  only  on  successive  measurements  of  position* 
The  state  transition  and  the  observation  matrices  are 

■  ■  [:;] 

H  =  Cl  03 

The  covariance  matrices  £*  P*  and  G  are*  of  course*  2  by 
2  matrices.  The  state  vector  has  two  elements*  velocity  and 
position*  while  the  measurement  vector  has  only  position* 
The  measurement  error  R  is  a  1  by  1  matrix  which  we  will 
vary*  We  hav<=:  chosen  G  as 


a  -  [s  ?] 

arguing  that  any  process  noise  will  be  contained  entirely  in 
velocity*  That  is*  there  can  be  no  random  motion  that  is  not 
caused  by  a  random  velocity*  Randomness  of  velocity  will  feed 
into  position  through  the  state  transition  matrix* 


Even  in  this  simple  case*  solving  analytically  for 
steady-state  Sain  in  terms  of  R  and  Q  requires  solving  a 
system  of  4th  order  polynomial  eauations.  We  will  opt 
instead  for  a  computer  solution.  The  reader  may  continue  to 
think  in  terms  of  the  noise  variance  ratio*  where  R  will  take 
on  the  values  1*  10*  100  and  Q  will  remain  constant  as  above. 
Since  there  is  only  one  non-2ero  term  in  the  Q  matrix*  we  may 
think  of  the  noise  variance  ratio  as  the  scalar  Quantity 
CK2*2)/R.  The  resulting  steady-state  gains  are 

Noise  Variance  Position  Velocity 


Ratio 

Gain 

Gain 

1.0 

*769 

♦  481 

0.1 

.553 

.211 

0.01 

.362 

.080 

As  would  be  expected*  the  position  gain  is  much  higher 
than  that  of  a  scalar  filter  at  an  eaui valent  noise  variance 
ratio*  because  the  process  variation  now  applies  to  velocity 
rather  than  position.  The  velocity  gain  is  considerably  less 

than  the  position  gain*  since  the  velocity  is  /not  measured 

/ 

t 

directly  but  must  oe  estimated  from  successive^  measurements 
of  position.  The  impulse- response  function  of  the  medium-gain 
filter  (noise  variance  ratio  0.1)  is  presented  in  figure  14. 

The  amplitude  response  of  these  three  filters  is 
compared  in  figures  15  and  16.  The  most  striking  feature  is 
the  amplification  which  occurs  at  a  specific  freeuency  in  the 
position  freouency  response.  This  implies  that  the  filter  is 
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VELOCITY  FILTER  IMPULSE— RESPONSE 
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FREQUENCY  RESPONSE  OF  POSITION  ESTIMATE 
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FRECtUENCY  RESPONSE  OF  VELOCITY  ESTIMATE 
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mast  sensitive  to  motion  in  a  particular  freouency  range* 
Thus *  the  natural  freouency  of  the  system  to  be  observed *  if 
it  is  known*  is  a  significant  design  parameter*  Again*  we 
observe  the  attenuation  of  high-f reouency  noise*  although  a 
significant  amount  still  remains  at  the  maximum  freouency 
(note  that  figure  16  does  not  include  the  origin)* 

The  freouency  response  of  velocity  shows  a  reduction  in 
amplitude  at  low  freouency.  The  amplitude  at  zero  freouency 
is  eoual  to  the  velocity  gain*  This  is  far  from  ideal 
performance  for  a  differentiator*  which  should  have  an 
amplitude  response  of  zero  at  zero  freouency?  with  a  slope  of 
one  up  to  the  cutoff  freouency  Cref*71.  The  differentiator 
is*  however*  reasonably  effective  at  reducing  the  amplitude 
of  high-f reouency  components. 

The  Phase  shift  of  the  filters  again  shows  increasing 
phase  lag  as  gain  is  decreased.  The  overall  effect  is 
similar  to  the  scalar  filter*  and  is  otherwise  unremarkable* 
Therefore*  plots  are  not  included. 

The  data  of  figure  8  was  tested  on  the  lowest-gain 
velocity  filter*  Note  that  the  spectral  content  of  the  data 
(figure  9)  is  ouite  iow-f reouencv;*  and  that  the  bandwidth  of 
the  lowest  gain  velocity  filter  is  ouite  wide*  and  indeed  is 
higher  than  that  of  our  scalar  filter*  So  it  might  be 
expected  that  the  velocity  filter  would  have  some  trouble 
with  the  data* 

The  velocity  filter  performance  on  the  nominal  function 
only  is  presented  in  figure  17.  The  filter  overshoots  ouite 
badly  at  the  discontinuous  steps*  which*  of  course*  are  an 
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impulse  in  velocity.  The  overshoot  is  less  at  the  start  and 
stop  of  the  trend.  It  does  settle  down  and  track  the  trend 
without  1 ad*  which  is  an  improvement  in  performance  over  the 
scalar  filter.  It  should  be  noted  that  with  hidher  dain*  the 
filter  would  track  the  nominal  data  better*  while  with  lower 
Sain*  the  overshoots  would  be  more  severe. 

The  performance  of  the  velocity  filter  on  the  <ata  is 
illustrated  in  fidure  18.  As  expected*  the  filter  tends  to 
follow  the  noise  too  much.  However*  it  does  follow  the 
discontinuities  much  more  Quickly  than  the  scalar  filter. 
This  points  out  the  fact  that  the  hidher-order  filter  is  more 
effective  as  a  maneuver  detector  but  it  is  less  suitable  for 
smoothind  very  noisy  data.  This  adain  illustrates  the 
concept  of  bandwidth*  which  is  auite  hidh  even  in  the 
low-dain  velocity  filter. 

In  retrospect*  the  decision  to  choose  8(1*1)  as  zero  may 
not  have  been  wise.  Allowind  some  process  noise  in  position* 
exclusive  of  velocity*  could  well  have  some  smoothind  effect 
on  the  velocity  estimate*  which  would  result  in  smoother 
one-period  ahead  predictions.  This  could  smooth  the 
operation  of  the  filter  a  bit.  The  possible  combinations  of 
filter  parameters*  even  for  this  simple  filter*  are  Quite 
numerous. 

The  freauency  response  of  a  second-difference 
(acceleration)  filter  was  also  determined  for  comparison. 
The  results  are  presented  in  fidure  19.  The  8  matrix  was 
zero  except  for  8(3*3)*  which  was  one.  R  was  chosen  as  10* 
resultind  in  3  nominal  noise  variance  ratio  of  0.1.  Adain* 
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VELOCITY  FILTER  PERFORMANCE 
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the  gain  was  higher  than  that  of  a  velocity  filter  with  an 
eouivalent  noise  variance  ratio*  The  amplif  icatiori  of 
low-f reouency  components  of  position  was  increased  and  the 
zero-f reouency  amplitudes  of  velocity  and  acceleration  again 
corresponded  to  filter  gain* 


FREQUENCY  RESPONSE  OF  ACCELERATION  FILTER 
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IV*  ESTIMATION*  SMOOTHING*  AND  PREDICTION 


U  ESTIMATION 


So  far*  we  have  been  concerned  only  with  estimation  of 
the  present  state.  A  filter  designed  to  provide  such  an 
estimate  cannot  be  symmetric*  because  it  can  put  no  weidht  on 
future  observations.  Thus*  phase  lad  is  inevitable*  and  is 
one  of  the  parameters  that  should  be  considered  in  the  desidn 
process. 


B.  SMOOTHING 


Smoothind  is  the  use  of  a  filter  to  provide  an  estimate 
of  past  states.  Such  a  filter  can  be  made  symmetric*  which 
completely  eliminates  the  Phase  lad.  Non-recursive  smoothind 
filters  cause  a  loss  of  N  data  points  at  each  end  of  the 
data*  where  the  span  of  the  filter  is  (-N*  N>» 

The  Kalman  filter  can  be  used  as  a  smoother  by  simply 
runnind  the  forward  estimate  throudh  the  filter  in  the 
opposite  direction.  The  impulse  response  function  of  the 
scalar  filter  was 


d<t>  =  ab 


It  can  be  shown  (appendix  A)  that  the  impulse  response 
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function  of  the  smoother  (forward  and  backward  filters 
combined)  is 


G(t)  =  ab*  /(1+b) 

which  is  Just  the  convolution 

g(t>  <D  g(-t) 

The  Kalman  filter  is  able  to  provide  an  estimate 
throughout  the  span  of  the  data*  No  data  is  lost  at  either 
end.  However*  due  to  transient  effects*  the  data  near  either 
end  is  subject  to  phase  shift  and  some  increase  in  gain.  The 
filter  is  necessarily  not  symmetric  near  each  end  of  the  data 
span. 

Gelb  Cref,43  includes  a  complete  discussion  of 
fixed-point*  fixed-lag*  and  fixed-interval  smoothing.  We 
will  restrict  our  attention  to  the  scalar*  fixed-interval* 
steady-state  case*  ignoring  the  end  effects. 

The  scalar  Kalman  filter  of  section  III.D  (noise 
variance  ratio  of  0.1*  gain  of  0.27)  was  used  as  a  smoother 
on  the  data  of  figure  8.  The  results  are  presented  in 
figure  20.  As  compared  to  the  one-pass  performance  as 
illustrated  in  figure  12*  the  smoothed  data  shows  phase  lag 
removed  and  peaks  in  the  oscillations  reduced.  However*  the 
smoother  has  less  ability  to  follow  the  discontinuities  in 
the  nominal  function.  The  removal  of  the  phase  lag  is 
characteristic  of  any  symmetric  filter.  However*  the  reduced 
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ability  to  follow  the  discontinuities  is  the  result  of 


reduced  effective  dein.  Since  we  have  convolved  the  filter 
weidhts*  we  have  squared  the  amplitude  of  the  freauenca 
response.  The  weidht  on  the  data  point  at  <t=0>  is  reduced. 
The  effective  dain  was  0.27  for  the  forward  filter*  and  0.156 
for  the  smoothind  filter.  This  reduction  in  effective  dain 
is  not  addressed  in  the  literature  on  the  Kalman  filter*  and 
it  is  unclear  how  smoother  dain  should  be  chosen  in  relation 
to  the  noise  variance  ratio  Q/R. 

C.  A  COMPARISON  OF  TWO  SMOOTHERS 

The  Kalman  smoother  of  the  previous  section  was  compared 
with  a  Gaussian  smoother  to  illustrate  some  desidn  options 
and  procedures.  The  Gaussian  smoother  was  chosen  from  amend 
a  hude  variety  of  data  windows  because  it  has  dood  smoothind 
properties*  and  because  it  is  particularly  easy  to  desidn. 
Interestindly*  preliminary  experiments  showed  repetitive 

i 

applications  of  a  Kalman  filter  to  result  in  an  approximately 
Gaussian  filter  weidht  distribution.  A  comparison  of  the 
Gaussian  smoother  to  a  variety  of  other  windows  is  contained 
in  Harris  Cref.153. 

The  Gaussian  smoother  is  a  symmetric  filter  with  the 
weidhts  chosen  accordind  to  a  discretized  and  truncated 
normal  distribution.  The  .formula  is 

d(t)  =  K  exp<-t*/2 C*) 


where  t  is  an  integer  on  the  range  (-*N*  N)  and  K  is  chosen 
such  that 

£  g(t>  =  l 

The  ease  of  design  comes  from  the  observation  that  the 
Fourier  transform  of  the  continuous  Normal  distribution  is 
also  a  Normal  distribution  with  scale  parameter  (variance) 
eoual  to  l/0*f  As  long  as  <<r  >2)  and  truncation  is  not  more 
severe  than  jNj  >  26>  a  reasonable  approximation  of  the 
freeuencu  response  for  the  Gaussian  digital  filter  is 

G(v)  Si  exp(-tf*v*/ 2) 

The  scale  parameter  was  chosen  such  that  the  freouenca 
response  was  eaual  at  <v  =  0*5) .  Skipping  the  algebraic 
details*  this  reouired  (&  -  3.11).  The  Gaussian  smoother  was 
truncated  to  <N  =  7>*  resulting  in  a  filter  span  of  15  data 
points.  The  freouenca  response  of  both  filters  is  presented 
in  figure  21*  and  the  filter  weighting  coefficients  are 
presented  in  figure  22. 

Since  we  truncated  the  Gaussian  smoother*  we  would 
expect  some  ripples  in  the  tail  of  the  freouenca  response* 
which  are  Just  barela  visible  in  figure  21.  The  Gaussian 
filter  has  a  sharper  transition  band*  and  is  Quite  effective 
in  blocking  high  freouencies.  As  compared  to  the  Kalman 
smoother*  the  Gaussian  filter  weights  the  present  data  point 
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FREQUENCY  RESPONSE  OF  SINGLE  AND  POUPLE  FILTERS 
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FIGURE  22 


less »  nearby  data  points  more  and  farther  data  points  less* 

The  performance  of  the  Gaussian  smoother  on  the  data  of 
figure  8  is  presented  in  figure  23*  A  comparison  of  the 
Kalman  and  Gaussian  smoothers  is  presented  in  figures  24  and 
25*  In  figure  24*  the  smoothers  are  applied  only  to  the 
nominal  function.  It  can  be  seen  that  the  Gaussian  smoother 
followed  the  discontinuities  and  corners  of  the  nominal 
function  better  than  the  Kalman  smoother*  However *  when  the 
smoothers  were  applied  to  the  noisy  data*  the  results  were 
less  clear  (figure  25).  The  Gaussian  smoother  again  followed 
the  nominal  function  a  bit  better*  but  it  also  followed 
low-f reouency  components  of  the  noise  a  bit  more*  tending  to 
emphasize  cyclic  effects  that  aren't  really  there.  The  mean 
souare  difference  between  the  smoothed  estimate  and  the 
nominal  function  were  very  similar*  1.81  for  the  Gaussian 
smoother  and  1.84  for  the  Kalman  smoother.  Thus*  the  Kalman 
filter  seems  ouite  effective  when  used  as  a  smoother.  The 
reader  is  reminded  that  the  mean  sauare  difference  for  the 
scalar  Kalman  filter  was  5.1*  which  clearly  indicates  the 
superiority  of  smoothing  over  filtering. 

D.  PREDICTION 

Prediction  is  difficult.  Recall  that  a  stochastic 
process  is  an  ensemble  of  possible  paths*  while  data  is  the 
manifestation  of  one  member  of  that  ensemble.  What  could 
have  happened  did*  but  what  can  happen  isn't  necessarily 
going  to.  Prediction  can  be  thought  of  as  filtering  without 
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GAUSSIAN  SMOOTHER  MKFORHAHCE 


COMPARISON  OF  GAUSSIAN  AND  KAL.HAN  SMOOTHERS 


measurements*  and  the  estimate  is  projected  forward  in  time 
through  the  state  transition  model*  Filtered  or  smoothed 
estimates  may  be  Quite  accurate*  even  if  the  state  transition 
model  is  not.  However*  good  predictions  are  heavily 
dependent  on  an  accurate  state  transition  model. 

The  higher-order  Kalman  filters  are  polynomial  models. 
Hamming  has  pointed  out  that  polynomial  models  are  poor 
predictors*  since  the  estimate  tends  to  veer  off  to  plus  or 
minus  infinity  as  soon  as  the  model  is  released  from  the  data 
Cref.73.  There  is  no  reauirement  to  use  the  same  model  for 
prediction  as  for  filtering.  For  example*  it  might  make  sense 
to  track  a  target  with  an  acceleration  filter*  but  to  compute 
fire  control  information  based  on  a  constant-velocity  model* 
since  target  acceleration  is  generally  assumed  random  with 
zero  mean.  Similarly*  the  economist  may  desire  to  filter 
data  with  a  high-order  model*  but  make  predictions  based  on 
constant  trend.  Clark  Cref.53  discusses  a  somewhat  more 
sophisticated  method  due  to  Singer*  in  which  the  model  decays 
from  an  acceleration  predictor  to  a  constant-velocity 
predictor  as  prediction  time  increases.  Such  techniaues  are 
heuristic  in  nature*  but  can  prove  valuable  to  the  innovative 
analyst. 

An  interesting  example  of  the  above  concept  can  be  found 
in  Box  and  Jenkins  Cref.lll.  They  compared  a  Quadratic 
forecast  due  to  Brown  Cref.103  with  their  own  IMA<0*1*1> 
model  with  a  gain  of  0.9.  The  latter  model  is  eauivalent  to 
the  steady-state  scalar  Kalman  filter.  Tk.»  data  used  for  the 
comparison  was  a  time  series  of  IBM  stock  prices.  Box  and 


Jenkins  observed  that*  while  the  Quadratic  model  roidht  well 
be  used  to  fit  the  data*  its  performance  as  a  predictor  was 
clearly  inferior  to  the  simpler  IMA(0*1*1)  model.  This  is 
not  surprising  si=v'e  it  has  land  been  suddested  that  stock 
prices  behave  as  a  random  walk*  and  that  the  best  forecast  of 
stock  price*  at  least  in  the  short  run*  is  the  present  price 
Cref.lll.  Note  that  the  foredoind  implies  that  the  dain 
should  be  set  to  1.0*  which  corresponds  to  no  filterind  at 
all.  Therefore*  Box  and  Jenkins  apparently  found  that  some 
filterind  of  the  data  was  appropriate*  even  thoudh  the  dain 

! 
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they  used  was  Quite  hidh. 
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V.  SOME  REFINEMENTS*  EXTENSIONS  AND  ALTERNATIVES 

A.  ADAPTIVE  FILTERING 

In  the  linear  Kalman  filter*  the  dai n  is  completely 
independent  of  the  data*  Clearly*  this  will  result  in  major 
errors  if  the  dain  is  chosen  inappropriately  or  if  the  data 
statistics  charide.  If  the  dain  is  too  low*  the  filter  lads 
badly*  In  the  extreme  case*  which  occurs  if  the  process 
noise  covariance  Q  is  much  too  low*  the  filter  pays  much  too 
much  attention  to  the  past  and  di verdes  from  the  data.  On 
the  other  hand*  if  the *dain  is  too  hidh*  the  filter  pays  too 
much  attention  to  the  data  and  the  state  estimate  contains 
noise.  If  the  filter  is  a  polynomial  model  and  is  to  be  used 
as  a  predictor*  the  resultind  errors  will  be  spectacular* 

The  solution  is  simple  in  concept  but  can  be  difficult 
to  implement.  One  simply  sets  the  steady-state  dain  as  low 
as  appropriate  for  the  stable  process  beind  estimated.  In 
tardet  trackind*  the  dain  would  be  set  to  track  an  airplane 
flyind  a  straidht  path.  A  'maneuver  detector*  or  'trend 
detector*  is  incorporated*  which  is  nothind  more  than  a 
recursive  statistical  test  applied  to  the  residuals  to 
determine  whether  or  not  they  come  from  a  zero-mean 
distribution.  If  not*  the  bandwidth  is  dradually  widened 
(dain  is  increased)  until  the  residuals  pass  the  zero-mean 
test.  Then*  the  dain  is  allowed  to  decrease  toward  the 
stable*  steady-state  value. 
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Further  details  and  some  novel  approaches  are  discussed 
in  Clark  E ref *53*  Two  examples  taken  from  Clark  are 
illustrated  in  fidure  26*  The  conceptual  adaptive  filter 
discussed  above  reouires  time  to  detect  the  maneuver  or 
trend*  adapt  to  it*  and  reconverde  to  a  stable  dain  settind* 
Durind  this  time*  th'  b^ate  estimate  is  less  accurater  and 
the  time  reetuired  maa  be  unacceptabla  lond  for  some 
applications* 

Clark  proposes  a  dual-bandwidth  adaptive  filter  to  speed 
adaptation*  The  process  is  simultaneousla  tracked  ba  a 
narrow-band  and  a  wide-band  filter*  If  a  maneuver  or  trend 
is  detected*  the  state  estimate  of  the  wide-band  filter  is 
fed  into  the  narrow-band  filter*  Ideally*  this  would  allow 
the  narrow-band  filter  to  Jump  immediatela  to  the  current 
(unbiased)  estimate  of  the  wide-band  filter*  In  practice* 
Clark  found  that  some  widenind  of  the  bandwidth  of  the 
narrow-band  filter  was  also  reoui red. 

Voluminous  literature  exists  on  the  subject*  much  of  it 
very  difficult  to  read.  Clark  Eref»53  incorporates  a 
particularly  lucid  account  of  stability  problems  encountered* 
methods  of  reducind  the  cost  of  false  detection  of  bias* 
analytical  methods  of  determiriind  filter  parameters*  and 
experimental  results.  Althoudh  Clark's  filter  was  desidned 
to  track  and  predict  the  position  of  airborne  tardets*  the 
methods  discussed  are  adaptable  to  the  filterind  of  economic 
time  series  or  virtually  any  other  stochastic  process. 
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B .  NON-LINEAR  FILTERING 


In  the  non-linear  Kalman  filter*  one  or  more  of  the 
matrices  Q»  R*  H*  or  <g  are  allowed  to  vary  with  time*  Since 
this  results  in  a  time-variation  of  the  sain  matrix* 
Quantitative  analysis  in  the  freouency  domain  is  no  lender 
appropriate*  However*  it  is  well  to  keep  the  concepts  in  mind 
in  order  to  gain  additional  insight*  There  are  two  basic 
types  of  non-linearities  that  may  arise*  non-linear 
measurements  and  non-linear  dynamics* 

1 •  Non-Linear  Measurements 

Non-linear  measurements  arise  when  observations  are 
made  in  one  coordinate  system  and  the  model  reouires  that  the 
state  be  estimated  in  another  coordinate  system.  In  this 
case*  the  matrices  R  and  H  are  time-varyins  functions  of  the 
coordinate  transformation*  and  do  depend  on  the  data*  in  the 
sense  that  they  depend  on  the  location  of  the  data  within  the 
coordinate  system.  This  type  of  non-linearity  is  often  easy 
to  handle. 

For  the  best  example  of  non-linear  measurements  we  must 
return  to  the  target-tracking  model*  Fire  control  systems 
generally  track  in  azimuth*  elevation*  and  range*  However* 
the  model  is  a  polynomial  in  Cartesian  coordinates*  but  not 
in  polar  coordinates*  Airplanes  often  fly  a  straight  path* 
but  seldom*  if  ever*  fly  a  constant  bearing  or  range  with 
respect  to  the  radar  observer*  In  this  model*  the 


non-linearity  can  be  reduced  by  considering  the  Cartesian 
measurement  error  as  a  linear  transformation  of  the  polar 
measurement  error*  If  the  polar  measurement  error  is 
Gaussian*  the  Cartesian  measurement  error  is  very  nearly 
Gaussian  with  covariance  matrix  R  a  function  of  the 
coordinate  transformation* 

The  Cartesian  R  matrix  will  not  be  diagonal*  even  if  the 
polar  R  matrix  is*  However*  Clark  Cref.53  has  found  that 
setting  the  of f-diasonal  terms  of  the  Cartesian  R  matrix  to 
zero  did  not  appreciably  degrade  filter  performance*  In  this 
way*  he  was  able  to  decouple  a  nine-state  filter  into  three 
three-state  filters. 

If  the  measurement  non-linearity  is  too  severe*  it  may 
not  be  reasonable  to  assume  that  the  noise  is  Gaussian. 
However*  limited  experiments  performed  on  data  with 
non-Gaussian  noise  (an  exponential  distribution  was  used) 
showed  that  the  Kalman  smoother  and  the  Gaussian  smoother 
were  Quite  robust  as  long  as  the  gain  was  not  high.  This 
seems  to  be  a  conseauence  of  the  Central  Limit  Theorem*  since 
low  gain  implies  a  linear  combination  of  a  fairly  large 
number  of  data  points.  It  should  be  noted  that  a  filter 
designed  to  handle  this  situation  is  still  linear*  although 
the  Gaussian  assumption  is  violated* 

2*  Non-Linear  Dynamics 

Non-linear  dynamics  are  considerably  harder  to  handle 
than  non-linear  measurements*  This  is  unfortunate*  since  the 
areas  of  potential  application  are  numerous*  Non-linear 
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dynamics  occur  when  the  Q  or  f  matrices  depend  on  the 
previous  history  of  the  process* 

As  a  simple  example*  consider  the  multiple  regression 

model 

y<t)  =  a  >^<t)  +  b  >^<t)  +  y<t-l) 

where  y<t-l)  corresponds  to  the  intercept  term*  In  this 
model*  we  wish  to  estimate  the  dependent  variable  y(t)*  To  do 
so*  we  need  to  estimate  not  only  the  independent  variables 
x/t)  and  }L<t>*  but  also  the  regression  coefficients  a  and  b* 
Let  us  assume  we  can  measure  y<t>»  x,(t)*  and  x^t )*  but  not  a 

and  b.  Assuming  a  first-order  system*  the  state  update 
eauation  is 
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Where  the  transition  matrix  is  unfortunately  not  unioue*  The 
first  row  of  the  transition  matrix  could  be  eeually  well 
represented  by 


C  1  0  0  Xj(  tit)  Xj<  t !  t )  3 
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C  1  a(t  i  t)/2  b ( t ' t >/2  x,<t!t>/2  xftlt)/2  1 

and  it  obviously  changes  at  every  iteration.  It  is  at  this 
point  that  filter  design  becomes  an  art. 

Note  that  the  independent  variables  and  the  regression 
coefficients  are  assumed  here  to  be  first-order  Gauss-Markov 
processes.  Increasingly  high  orders  would  multiply  the  state 
space. 

Several  experiments  were  run  using  a  second-order  model 
similar  to  the  above  on  the  Box- Jenkins  C ref. Ill  series  M 
data  (sales  data  with  leading  indicator).  Quantitative 
results  are  not  presented*  because  the  Box-Jenkins  data  did 
not  include  sufficient  forecast  estimates  for  comparison* 
some  "cheating*  was  done  because  the  Box-Jenkins  parameters 
were  used  in  filter  design*  and  it  never  became  clear  exactly 
what  parameters  were  appropriate  for  the  R  and  Q  matrices. 
However*  some  Qualitative  comments  are  appropriate.  The 
model  did  work.  Some  instability  was  noted  in  the  regression 
parameters.  It  became  obvious  that  the  gain  on  the 
regression  parameters  must  be  set  very  low  in  comparison  to 
the  gain  on  the  independent  variables*  in  order  to  keep  the 
regression  parameter  estimates  from  varying  faster  than  the 
estimates  of  the  independent  variables.  This  implies 
choosing  small  values  for  the  noise  variance  of  the 
regression  parameters.  Also*  by  keeping  the  gain  fairly  low 


on  the  leading  indicator*  it  was  possible  to  induce  a  phase 
lad  that  approximately  cancelled  the  lead* 

The  reauirement  to  keep  dain  low  in  order  to  improve 
stability  is  evidently  a  conseauence  of  the  increased  decrees 
of  freedom*  The  more  parameters  to  be  estimated*  the  more 
dedrees  of  freedom  in  the  model*  Hidh  dain  is  analadous  to 
relatively  few  data  points  beind  used  in  redression.  The 
more  variables  we  introduce  into  the  model*  the  less  dain  we 
are  able  to  use* 

It  is  indeed  unfortunate  that  multiple  redression  is  a 
non-linear  problem  when  cast  in  a  filter  model*  It  would  be 
useful  +o  have  a  multiple  redression  model  for  which  more 
recent  observations  were  weidhted  more  heavily  than  older 
ones  in  determinind  the  redression  parameters.  No  doubt  the 
innovative  analyst  could  develop  one  to  fit  the  specific 
situation.  However*  clearly-defined  technioues  with 
demonstrated  results  are  not  available  to  the  practitioner* 
The  experts  all  have  their  favorite  methods*  and  much  of  the 
literature  is  difficult  to  read*  There  is  clearly  a  need  for 
additional  research  in  this  area. 

C.  NQN-PARAMETRIC  FILTERING 

We  close  our  discussion  with  an  interestind  alternative 
to  conventional  didital  filterind  technieues.  There  are 
those  who  are  bothered  by  the  usual  distributional 
assumptions  made  in  any  application  of  parametric  statistics. 
An  extensive  literature  has  developed  in  the  field  of 
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non-parametric  statistics*  which  is  based  on  the  principle 


that  distributional  assumptions  are  avoided*  or  at  least 
weakened.  A  strand  paint  of  non-parametric  statistics  is 
relative  insensitivity  to  extreme  outliers.  However*  little 
prodress  has  been  made  in  the  non-parametric  analysis  of  time 
series.  An  exception  may  be  found  in  Tukey  Cref.163*  which 
is  presented  in  a  hidhly  intuitive  manner*  with  little  or  no 
theoretical  backdround. 

One  simple  idea  advanced  by  Tukey  is  that  of  median 
smoothind.  The  smoothed  estimate  i  ■>  based  on  the  median  of 
several  adjacent  data  points*  rather  than  on  a  weidhted 
linear  combination.  The  result  is  ouviously  a  series  of 
steps*  since  adjacent  data  points  will  often  have  the  same 
median.  Tukey  susdests  several  methods  to  restore  some 
curvature  in  the  estimate.  The*  a  will  not  be  developed  here. 
Tukey 's  methods  would  be  relatively  hard  to  mechanise  on  a 
computer*  because  the  methodolody  reouires  extensive  lodical 
rules. 

Tukey ' s  methods  could  be  extended  to  real-time  filterind 
problems  by  developind  a  non-parametric  analod  to  the 
recursive  dxdital  filter.  Recall  that  the  recursive  didital 
filter  consists  of  a  weidhted  linear  combination  of  recent 
data  points  added  to  a  weidhted  linear  combination  of  recent 
estimates.  Th  •  non-parame<vric  filter  estimate  would  simply 
be  the  median  of  several  recent  data  points  and  several 
recent  states.  The  idea  is  intuitively  appealind  and  should 
be  the  subject  of  future  research.  Discussion  here  will  be 
limited  to  some  of  the  more  obvious  traps  that  await  the 
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unwary 


If  such  a  filter  were  to  be  designed*  the 
impulse-response  function  would  he  meaningless *  because  the 
median  estimate  would  always  be  zero  if  the  span  of  the 
filter  were  greater  than  two*  Nevertheless*  a  median  filter 
does  have  a  freauency  response*  which  in  fact  is  a 
particularly  nasty  one* 

Consider  a  seven-point  median  smoother*  where  the  state 
estimate  at  time  t  is  the  median  of  the  measurements  made  at 
time  <t-3)  to  time  (t+3>*  This  is  analogous  to  the 
rectangular  (parametric)  window  discussed  in  Hamming  C ref. 73* 
The  rectangular  window  weights  all  data  points  within  the 
window  eoually*  The  median  window  obviously  does  the  same* 
As  a  result*  we  would  expect  the  freauency  response  of  the 
median  window  to  have  severe  ripples  as  does  that  of  the 
rectangular  window*  We  can  see  intuitively  that  this  is  true* 
Since  the  span  of  our  example  median  window  is  7*  the 
freauency  response  of  any  freauency  that  is  a  non-zero 
integer  multiple  of  1/7  is  obviously  zero.  The  freauency 
response  at  zero  freauency  is  one*  since  the  zero  freauency 
implies  a  constant*  The  amplitude  of  the  freauency  response 
falls  off  to  the  first  zero*  then  rises  again*  Successive 
maxima  decrease  with  increasing  freauency*  but  the  freauency 
response  is  always  non-zero  except  at  freauencies  that  are 
non-zero  integer  multiples  of  the  reciprocal  of  the  span  of 
the  window*  Thus*  the  non-parametric  filter  will  need  to 
incorporate  some  (unspecified)  device  to  improve  the 
freauency  response* 
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Two  other  difficulties  are  worthy  of  mention*  First* 
the  sampling  distribution  of  the  median  may  have  a  larder 
variance  than  the  sampling  distribution  of  the  mean*  This 
means  that  the  parametric  filter  may  provide  a  better 
estimate  than  the  non-parametric  filter  if  the  assumptions  on 
which  the  parametric  r.J.ter  is  ba-^d  a'  'v  all  reasonable* 
Second*  for  the  n;  --*•  ^  filter  o  be  useful*  the 

median  must  be  a  of  interest*  If  it  is  assumed 

that  the  distribu-  etric*  the  median  and  r>ean  3re* 

of  course*  eoual*  .  .  sam.-1’.  rg  distribution  is  sK  : 

mean  cannot  be  cleducec,  from  the  median  unle  rict 

parametric  assumptioris  are  imposed*  which  of  course*  override 
the  Justification  for  the  non-pararoetric  filter  in  the  first 
place.  The  idea  is  nevertheless  intriguing*  3nd  should  be 


explored  further 


VI.  SOME  APPLICATIONS 


The  Kalman  filter  has  been  applied  to  Operations  Research 
and  economic  problems  with  varying  degrees  of  success. 
McWhorter  Href .173  conducted  an  empirical  study  of  the  Kalman 
filter  in  which  he  compared  it  to  several  other  methods  of 
time  series  forecasting.  The  results  were  mixed;  with  no 
method  dominating.  The  Kalman  filter  compared  more  favorably 
over  a  short  term  forecasting  horizon  than  over  a  long  term 
one.  Its  performance  was;  not  surprisingly;  found  to  be 
degraded  if  the  structural  model  was  seriously  mis-specif ied. 
McWhorter  pointed  out  some  of  the  difficulties  encountered  in 
building  the  model.  In  an  economic  context;  it  is  often  very 
difficult  to  specify  the  noise  covariance  matrices  R  and  Q; 
and  even  to  identify  the  structure  of  the  state  transition 
matrix  The  assumptions  made  are  often  sweeping  and 
arbitrary;  in  contrast  to  tracking  applications  where  the 
noise  processes  and  especially  the  state  transition  model  are 
relatively  well  understood. 

A.  INVENTORY  MANAGEMENT 

The  Kalman  filter  is  directly  applicable  to  inventory 
management;  and  if  properly  designed;  should  be  superior  bo 
the  finite  exponential  smoothing  model  of  Bessler  and  Zehna 


Cref.143.  Downing;  Pike;  and  Morrison  Cref.183  designed  a 


Kalman  filter  for  the  inventory  control  of  nuclear  material. 
The  paper  is  readable*  and  the  filter  is  well-documented  and 
easy  to  understand.  They  use  the  concept  of  a  control 
vector*  which  has  not  been  mentioned  here.  An  interesting 
peculiarity  of  the  model  is  that  one  of  the  measurements  is 
only  available  once  every  twenty  iterations.  The  state 
transition  matrix  is  a  simple  material  balance  relation  which 
is  obviously  Quite  accurate.  Such  a  model  could  be  expected 
to  perform  Quite  well. 

B.  ESTIMATING  A  MEAN  FUNCTION 

i 

Although  the  Kalman  filter  was  derived  from  an  assumption 
of  stationer ity *  we  have  seen  that  it  can  be  Quite  powerful 
in  separating  a  time  varying  signal  from  noise.  The  examples 
of  section  III  were  all  essentially  estimates  of  the 
time-varying  mean  function  of  a  stochastic  process.  The 
example  process  was  Gaussian  with  a  constant  variance.  The 
variance  was  the  measurement  noise*  and  so  directly 
influenced  the  gain*  If  variance  were  not  constant*  the 
performance  of  a  nan-adaptive  filter  would  be  degraded.  If 
the  change  in  variance  was  great  enough*  an  adaptive  filter 
would  be  reauired. 

A  good  method  of  estimating  a  time-varying  mean  function 
could  be  applied  in  numerous  areas *  such  as  any  sort  of 
traffic  or  flow  control  problem*  perhaps  in  Quality  control 
of  large-batch,  or  flow  manufacturing  processes*  and  any 
application  where  it  is  desirable  to  detect  a  change  in  the 


process*  The  sensitivity  of  the  filter  is  directly 
adjustable  by  the  modeller  through  the  noise  covariance 
matrices  Q  and  R. 

A  particularly  useful  application  would  be  to  the 
estimation  of  the  rate  parameter  of  a  non-time-homogeneous 
Poisson  process*  If  this  can  be  done  accurately*  the  process 
can  be  transformed  to  a  stationary  one  !Iref*63*  which  greatly 
expands  the  number  of  analytic  tools  that  can  be  used* 

The  Poisson  process  is  a  counting  process  in  continuous 
time*  and  to  attempt  to  filter  a  string  of  interarrival  time 
data  would  violate  the  sampling  theorem*  The  times  of 
arrival  are  the  measurement  times*  and  they  are  most 
certainly  not  made  at  eeually  spaced  intervals.  Instead*  the 
filter  may  be  designed  to  sample  a  counting  process.  At 
discrete  intervals  the  filter  would  count  the  number  of 
arrivals  since  some  arbitrary  time  origin.  If  the  process 
were  to  continue  for  a  long  time*  the  time  origin  might 
occasionally  have  to  be  reset  to  prevent  computer  overflow. 
It  is  eas  •  to  see  how  this  sampling  process  could  be 
implemented  even  if  the  input  data  were  actually  arrival 
instants  in  continuous  time.  The  sampling  interval  should  be 
small  enough  that  there  is  low  probability  that  more  than  one 
arrival  would  occur  during  a  given  measurement  interval. 
Since  the  number  of  arrivals  is  monotone  non-decreasing  in 
time*  a  velocity  or  trend  model  would  be  appropriate.  The 
input  data  would  consist  of  integers.  The  state  estimates 
would  not.  The  non-integer  estimate  of  number  of  arrivals  up 
to  the  current  time  would  not  be  useful  to  us.  However*  the 
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second  element  of  the  state  vector*  the  velocity  or  trend* 


would  in  fact  be  the  filtered  arrival  rate  estimate.  Since 
the  process  is  noisy  and  non-Gaussian*  a  very  low 
steady-state  dain  is  appropriate. 

The  time-varyind  Poisson  process  cannot  have  constant 
variance*  since  the  mean  and  the  variance  are  equal.  A  low 
arrival  rate  implies  hidh  variance  in  the  Poisson  process* 
which  is  equivalent  to  hidh  measurement  noise*  which  requires 
low  dain.  A  constant-dain  filter  would  therefore  be 
relatively  more  sensitive  at  low  arrival  rates  than  at  hidh 
arrival  rates.  An  adaptive  filter  could  be  easily  desidned 
to  use  the  inverse  of  the  rate  estimate  as  the  measurement 
noise  variance  estimate.  Stability  midht  require  that  the 
adjustment  of  the  measurement  noise  variance  be  itself  a 
filterind  process*  in  which  the  incomind  variance  estimate  is 
redarded  as  data. 

C.  MULTIPLE  REGRESSION 

If  the  redression  constants  are  assumed  known  (or 
computed  by  other  means)  the  desidn  of  an  appropriate  filter 
is  quite  straidhtforward*  and  quality  of  estimation  is 
related  directly  to  the  quality  of  the  model.  Note  that  the 
velocity  filter  is  simply  the  redression  of  velocity  on 
position*  where  the  slope  parameter  is  known  to  be  one.  If 
the  redression  coefficients  are  assumed  to  vary  in  time*  the 
problem  becomes  non-linear  and  is  quite  complex.  Because  of 
the  immense  applicability  of  this  model*  additional 
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developmental  work  is  indeed  a  fertile  field  for  future 
research* 

D*  SOME  DESIGN  CONSIDERATIONS 

In  applications  where  the  noise  covariance  matrices  R  and 
0/  and  the  system  dynamic  model  (state  transition  matrix  f) 
are  known  or  easily  estimated/  design  is  straightforward 
and  has  been  successfully  accomplished  while  remaining  in  the 
time  domain*  However/  in  applications  where  sweeping 
assumptions  are  reauired/  a  f reauency-domain  analysis  could 
be  very  helpful.  Some  guidelines  are  as  follows! 

1 .  Spectral  Analysis  of  the  Data 

A  spectral  analysis  of  sample  data  will  s'  ou  what  the 
freauency  response  of  the  filter  should  be.  The  Fa<=:  Fourier 
Transform  (,FFT)  program  available  in  most  computer  libraries 
is  generally  easy  to  use.  However/  the  FFT  programs 

J 

generally  reeuire  an  exact  power  of  2  for  the  number  of  data 
points.  Hamming  Cref.73  points  out  some  pitfalls.  Since 
stationarity  is  assumed/  the  data  should  be  considered  as  a 
rotating  cylinder/  and  if  the  starting  and  ending  values  are 
not  similar/  a  discontinuity  will  exist  in  the  spectrum.  The 
data  can  be  tapered  and  padded  with  zeros/  but  exactly  the 
best  method  to  accomplish  this  is  unknown*  Several  methods 
might  be  tried. 

The  main  virtue  of  the  FFT  is  its  speed.  It  works  well 
on  a  long  run  of  data.  If  the  number  of  data  points  is  small 
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(around  a  hundred)  it  might  be  effective  to  find  (or  write)  a 
less  efficient *  conventional  discrete  Fourier  transform 
program*  which  would  not  require  padded*  truncated*  or 
tapered  data  if  the  starting  and  ending  values  are  similar* 
If  the  FFT  program  used  does  not  reauire  an  exact  power  of 
two  for  the  number  of  input  data  points*  it  would  be  well  to 
find  out  why  not*  The  program  may  be  doing  the  padding  and 
tapering  itself*  and  the  analyst  should  be  curious  as  to  how. 
The  analyst  should  remember  that  the  spectrum  is  computed 
from  the  data*  and  it  is  therefore  an  estimate.  If  the  run 
of  data  is  short*  there  will  be  considerable  variance  in  the 
estimate* 

2*  Frequency  Analysis  of  Proposed  Models 

The  analyst  may  test  the  effect  of  assumptions  made 
in  designing  the  filter  by  simply  obtaining  an  impulse 
response  of  the  filter  and  running  it  through  an  FFT. 
Truncation  and  tapering  is  no  problem*  because  the  impulse 
response  will  approach  zero  with  time.  The  proper  impulse 
function  is  simply  a  1  followed  by  2n-l  zeros  for  a  filter* 
or  a  1  in  the  middle  of  2*  -1  zeros  for  a  smoother.  If  the 
output  of  the  FFT  consists  of  real  and  imaginary  components* 
it  will  be  necessary  to  compute  amplitude  and  phase. 

3.  Adjusting  the  Model 

If  the  model  dynamics  seem  adequate  but  the  bandwidth 
is  wrong*  the  analyst  should  by  now  have  some  insight  into 
what  adjustments  to  make  to  the  noise  covariance  matrices  to 


try  to  improve  things.  In  3  model  of  any  complexity  at  all* 
there  are  numerous  possible  combinations*  However*  even  some 
improvement  over  the  initial  assumptions  will  be  beneficial* 
We  are  not  looking  for  theoretical  elegance*  we  are  looking 
for  performance. 

Perhaps  the  model  dynamics  obviously  call  for  a  trend 
filter  or  even  a  chande-of -trend  (acceleration)  filter*  but 
the  data  is  au-te  noisy.  Consideration  should  be  given  to 
lowering  the  order  of  the  filter*  A  very  low-gain  velocity 
filter  will  not  follow  changes  in  trend  well.  A  higher-gain 
scalar  filter  may  do  so  more  effectively*  although  it  will 
lag  a  steady  trend*  There  are  many  tradeoffs*  and  we  cannoT. 
achieve  perfection* 

4.  Testing  the  Model 

The  model  should  be  tested  on  real  or  simulated  data* 
From  here  on*  the  modelling  process  is  the  standard  cyclical 
one*  going  back  to  earlier  steps  as  necessary  until 
satisfactory  performance  is  achieved. 
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APPENDIX  A.  DERIVATIONS 


1 ♦  SCALAR  KALMAN  FILTER 


3*  Recursive  Formula  for  Kalman  Gain 

The  covariance  extrapolation  eouation 
P<t>  =  J  £Ct-i>  jgT  +  Q 
reduces  in  the  scalar  case  to 

P(t>  -  2<t-l)  +  Q 

Since 

K(t>  =  2(t>  htr*' 

we  may  writer  for  the  scalar  case* 

2<t-l)  *  K<t-1 >  R 

Similarly*  since 

t  r  *1 

K< t )  =  PH  C  HPH  +  RD 

by  reducing  to  the  scalar  case  and  substituting*  we  may  write 

K< t )  =  K<  t-1 )  R  +  Q  »  K(t-l)  +  Q/R 

K < t-1 )  R  +  Q  +  R  K<t-l )  +  Q/R  +  1 


b.  Steady-State  Kalman  Gain 

Rearranging  the  recursive  slain  eouation  and  letting 
K<t)  =  K'(t-l)  =  K 

we  see  that 

K*+  <Q/R)K  -  Q/R  *  0 
By  the  Quadratic  formula* 


I  X 

Q 

J  Q  + 

/  4R* 

R 

We  are  obliged  to  take  the  larger  root*  since  the  smaller 


root  would  force  the  dain  to  be  negative.  We  also  observe  the 
inverse  relationship 


<K  +  Q/2R>*  =  Q*/4R*+  Q/R 
kS  <Q/R>K  +  Q1/4R1  =  Q*/4R*“+  Q/R 
K*  =  <Q/R> < 1~K) 


Q/R  =  K  /<!-«> 


c.  Transient  Kalman  Gain 

Recall  that  the  Kalman  filter  reouires  a  prior  state 
estimate  X(0)  and  a  prior  estimate  of  covariance  P<0).  This 
reauirement  can  be  avoided  by  usind  K(l>  =  If  which  allows 
the  initial  state  estimate  to  be  eaual  to  the  first 
measurement.  Recall  that 

K< t)  =  I<t>  HTR*J 


Since  K<1)  *  1*  then  2<1>  =  R 


d.  Amplitude  and  Phase  of  Freauoncy  Response 


The  freauency  response  is 


since 


then 


H<v)  =  a  ]^Cb  exp<-iv> 
tso 


}b  e;,'P<-iv>|  <  1 


H<v)  =  a  /  Cl-b  ewp(~iv)3 


The  amplitude  souared  is 

A*  =  H<v)  H<~v>  =  a4/  Cl-b  e;<p(-iv)3  Cl-b  exp<iv)3 
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By  Euler's  relation 

A*  «  a*/  Cl  +  b*  -2b  ccs  vl 
which  may  be  written 

A*  a  3*/  C(l-b)*  +  2b (i -cos  v)3 

recalling  that 


Q/R  «  K* /( 1-K)  =  a2/b 
we  may  write  the  amplitude  as 


A 


Q/R _ 

Q/R  +  2 (1-cos  v) 


The  phase  andle  is 

0(v)  =  arctan  CIm(v>  /  Re(v)l 

whpfe  Im(v)  and  Re(v)  are  the  imaginary  and  real  par ti  of 

H(v)»  which  may  be  written 

H(v)  =  a _ Cl-b  exp(iv)3 

Cl-b  exp(-iv)  HCl-b  exp(iv):! 

H(v)  =  a  ( 1  frb  cos  v  -ib  sin  v) 

1  f  b*  -2b  cos  v 

which  allows  us  to  write 

0(v>  =  arctan  C(-b  sin  v)  /  (1-b  cos  v)3 

The  aridle  for  maximum  phase  shift  occurs  when 

d0  ( v)  =  fa*  -  b  cos  v  a  o 

dv  1  +  bx  -2b  cos  v 

so  that  the  maximum  phase  shift  0(v)max  occurs  when 

v  =  arccos  b 

and  has  a  value  of 

0<v)m3x  a  arctan  -fa  sin(arccos  b) 

1-b  cos (arccos  b) 

=  arctan  (  -b  /  J  1-b*  > 
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2, 


IMPULSE  RESPONSE  FUNCTIONS 

The  impulse-response  of  the  scaler  Kalman  filter  is 
a<t)  =  ab*  t  =  Or  l?2r.,. 

a*  Impulse-response  of  Double  Filter 

si<t>  =  ab  ®  ab 

40  * 

. . .  .  <  J*  Jt~k  *  ^  .4  . 

S(t)  =  £  ab  ab  =  a  ^£,b  b 

2<t)  =  (t+1)  a,b*‘ 

b.  Impulse-response  of  Scalar  Kalman  Smoother 

s(t)  =  ab*  ®  ab  * 


sKt)  =  ^.a  ab**^  =  a2j£b**fc=  a*b* 

A*-» 


a<t)  =  a*b*/  <l-b*>  =  a  b*/a+b> 
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